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List of Symbols 


The following is a list of all symbols used in this study, 
accompanied by a brief description. Vectors are indicated by 
brackets with the dimensions of the variable in the brackets. 
Matrices are denoted by brackets with both dimensions given 
within the brackets. 
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I . Introduction 


The design procedure for any complex system is typically 
one of compromise and trade-off. The design of the space 
shuttle reentry vehicle is certainly no exception because 
numerous factors demand attention and contribute to the per- 
formance capability of the vehicle. 

A primary consideration in the design of any reentry 
vehicle is the aerodynamic heating which the vehicle will en- 
counter upon reentry. The total heat input and the associated 
temperatures directly determine the amount of thermal protec- 
tion necessary for the safe reentry of the vehicle. The amount 
of thermal protection required is particularly important because 
the resultant weight is usually a significant portion of the 
vehicle weight. Consequently, weight penalties incurred by ex- 
cessive heating reduce the payload capability with a resultant 
increase in payload delivery cost. Therefore, trajectories 
which yield minimal weight penalties due to heating effects are 
desirable in order to improve the operational efficiency of the 
vehicle . 

Of the suggested models for the Thermal Protection System 
(TPS) , the Reuseable Exterior Insulation system (REI) is the 
model selected for this study. The REI system is composed of 
titanium covered by a surface insulation material. The insula- 
tion is used on that portion of the surface where the local 
temperature exceeds 650° F. Additional internal insulation is 
utilized on all fuselage and wing areas which exceed 650° F. 

To estimate the weight of the TPS, the surface of the ve- 
hicle is divided into n panels. The weight of each panel is 
related to the total heat transferred to the panel during re- 
entry. A typical variation of exterior insulation weight per 
unit surface area with heat load is indicated in Figure 1. 
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Heat Load/Unit Area (BTU/FT 2 ) x 10“ 3 


Figure 1. REI Weight Versus Heat Load 
If the curves shown above are approximated by a linear function 


as indicated by the dashed line, the TPS weight can be related 


to the heat 

load to the 

n panels as follows: 




n 



W TPS 

E w. 
i-1 1 

(1.1) 

where 

w . = 
1 

A + B Q- 

W W X 1 

(1.2) 



n 


or 

W TPS 

E (A + B Q.) 
w w x i J 

(1.3) 

A and B 

w w 

are constants 

while is the heat load to the 



i-th panel. Consequently, with these simplifications the TPS 


weight is linearly related to the heat load, and the minimiza- 
tion of the TPS weight is analogous to minimization of the total 
heat load. 
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II. Performance Index 

In the present study only the four lower-surface panels of 
the fuselage are considered. In addition, for reentry veloci- 
ties on the order of earth orbital speed, aerodynamic heating 
is composed primarily of convective heating. 1 Therefore, the 
heating rate to the underside of the fuselage is represented by 
the sum of the convective heating rates to the four individual 
panels, and the performance index to be minimized is the inte- 
grated heat input to these panels. 


J = 


X 


c 



( 2 . 1 ) 


where A is a constant scale factor and 6 is the convective 
c x c 

heating rate to the selected panels, given by 


Q 


c 


4 

q £ s . y . g . 

° i=i 1 1 1 


( 2 . 2 ) 


The reference heating rate, q Q , is the heating rate which 
would occur at the stagnation point on a one-foot radius sphere 
following the same path as the vehicle, that is, 


= k p°- 5 v 3 * 15 
c 


( 2 . 3 ) 


where k c is a constant, p is the atmospheric density and V 
is the magnitude of the velocity. The area of the i-th panel is 
designated by s^ . The function g^ is a boundary layer de- 
pendent function which relates the effects of the nature of the 
flow in the boundary layer on the heating rate to the i-th panel. 
It is defined^ to be the following second-derivative-continuous 
function : 
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g; " 1, 
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R 


L. 
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< 1.0 (Laminar) 


g- 


*i 


■f.o Vi 


R 

l 


R 

J 1 0 < - 2 - 
’ 1 * u - R t 


£ 1.5 (Transition) (2.4) 


R 


= k 


g III. * V R 


C^-) 


0.3 


R 


n 


r l. - 
1 


> 1.5 (Turbulent) 


n 


where R_ is the reference Reynolds number per foot 

r2 


R = k 
n 


R 


(7 


pV‘ 


-) 


1.5 


q o !/2 


In these formulas 


"R’ 


a . 
J 


are constants and R. 


(2.5) 


is the 


boundary layer transition Reynolds number for the i-th panel. 
The function y^ is introduced to account for the effect of 
angle of attack on the heating rate to the i-th panel. The par 
ticular form of y^ was chosen because it expresses the ex- 
pected influence of a and was simple to implement. The ex- 
pression for y. is given as 


y . (a) = b . + c - sin' 

1 1 y J l i 


a 


( 2 . 6 ) 


where 


and 


c^ are constants determined by fitting y. to 


experimental data and 


a 


is the angle of attack. 


III. Dynamical Model 

The dynamical model for the atmospheric reentry is chosen 
to be as uncomplicated as practical and yet retain the salient 
characteristics of a more exact formulation. In particular, the 
earth is assumed to be spherical, non-rotating and to possess an 

inverse square gravitational force field. The atmosphere is 
considered to be at rest with respect to the earth and to vary 
exponentially with altitude. 
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The spacecraft is located relative to the earth via a 
spherical coordinate system whose origin is fixed at the cen- 
ter of the earth as illustrated in Figure 2. The distance from 
the center of the earth to the vehicle is designated by r , 
while the longitude and latitude of the vehicle are represented 
by 0 and <f> , respectively. 

The velocity vector, with magnitude V , is oriented in 
space by using the heading angle, ip , and the flight path 
angle, y as indicated in Figure 3. The attitude of the ve- 
hicle is then established by a roll angle, p , about the veloc- 
ity vector followed by an angle of attack, a , about the ve- 
hicle's lateral axis as shown in Figure 4. Zero sideslip angle 
is assumed. 

In this system, the equations governing the motion of the 
vehicle are^ 


r = 

0 = 

<i> = 

v = 


V sin y 

V cos y cos \p 
r cos <j> 

V cos y sin }p 

r 

k sin y _ D 


Y = 

i> = 


k cos y # V cos y 
Vr 2 r 


V cos y cos ip tan 4> 
r ~ 



L sin p 
mV cos y 


(3.1) 


where k is the gravitational constant of the earth, m is the 
mass of the vehicle, L is the lift of the vehicle, D is the 
drag of the vehicle and motion of the vehicle about its center 
of mass is ignored. Lift and drag are defined in the conven- 
tional manner as 

L = (1/2) pV 2 SC L 
D = (1/2) pV 2 SC D 


and 


(3.2) 
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Figure 4. Definition of Control Angles 
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where S is the; vehicle reference area and Cj and C () are 
coefficients o f lift and drag, respectively. The atmospheric 
density, p , is assumed to vary exponentially with altitude 
according to 

P - P 0 - k d( r ' r e) (3.3) 

where p Q and are constants chosen to approximate the 

density of the actual atmosphere over the altitudes of interest 

for reentry and where r g is the radius of the earth. 

The aerodynamic coefficients and are generally 

functions of the Mach number and of the Reynolds number as well 

as the angle of attack. However, for hypersonic flight, the 

drag coefficient is essentially independent of the Mach number 

and for high altitude flight the effects of Reynolds number are 

relatively unimportant in comparison with those due to angle of 

attack. Consequently, the lift and drag coefficients are 

assumed to be functions of angle of attack only and are given by 

8 9 

the following relationships ’ obtained from Newtonian flow 
theory : 


sin a cos a |sin a| 
o 

(3.4) 

= + 1 s ^ n a l 

o L 

(3.5) 


C„ and Cp are constants. The lift-to-drag 
o U L 

is then given as 

sin a cos a |sin a| 

E = — 5 = (3.6) 

C D + Cp sin -3 a| 
o U L 

The dependence of C L , and E upon angle of attack for the 

vehicle of this study is illustrated in Figure 5. 


where , 

o 

ratio, E , 
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Figure 5. Aerodynamic Coefficients versus Angle of Attack 











The reentry vehicle is controlled by varying the angle of 
attack, a , which determines the magnitude of the aerodynamic 
force and the roll angle, 3 , which determines the direction 
of the lift force. Although for the analysis of this study, f? 
is not necessarily subject to physical limitation, a must 
certainly be limited. The desire to reduce heating on the upper 
surface of the vehicle requires that a be non-negative. In 
addition, the obvious constraint that a not exceed 90 degrees 
must be imposed. However, more exact modeling requires even 
closer restrictions to be placed on the angle of attack. For 

l 

example, control power and stability boundaries for a typical 
configuration place additional constraints upon the angle of 
attack as illustrated in Figure 6. 
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Figure 6. Angle of Attack Boundaries 
versus Mach Number 
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In accordance with such restrictions, the angle of attack 
limits are chosen to be 


a v = 55° 
max 


a . =15° 

min 


(4.1) 


Note that the maximum angle of attack is approximately equal to 
the angle of attack for maximum lift coefficient. 

Minimization of heat load is usually accomplished by pro- 
ducing high peak heating rates over short time intervals. This 
is done by flying a trajectory which is composed of a sequence 
of skip segments in which the vehicle dives down into the atmos- 
phere until sufficient lift is generated to force it back up into 
the thinner atmosphere. For a vehicle constrained to positive 
angles of attack, the skip is produced by first rolling down 
( | B| _> 90°) . This produces a downward component of lift which 
forces the vehicle down into the atmosphere. The pull-up is then 
accomplished by rolling up. Skipping maneuvers such as these can 
produce high peak altitudes as well as high heating rates requir- 
ing increased accuracy from the guidance system. As a means of 
reducing the adverse effects of the skipping phenomenon, a maxi- 
mum roll angle is implemented to reduce the vehicle's ability 
to dive into the atmosphere. Therefore, the constraints applied 
to the roll angle are 


e 

max 

^min 


70° 

0 ° 


(4.2) 


V, Boundary Conditions 

The immediate projected use of the shuttle vehicle is that 
of transporting men and equipment to orbiting earth satellites 
or space laboratories. Therefore, the initial conditions for 
reentry will remain approximately the same for most of its 
missions. The nominal values selected for this investigation 
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are indicated below: 

r - r = 400,000 ft 
o e 

e o = o 

*0 = 0 

V Q = 25,975 ft/sec (5.1) 

Y 0 = - 1.5 deg 

4> 0 = 0 

*o = 0 

The shuttle vehicle represents, in many respects, a signi- 
ficant departure from previous operational reentry vehicles. 

The increased hypersonic lift-to-drag ratio (1.5 to 2.0) of the 
shuttle enables it to achieve a larger reentry footprint. In - 
addition, the utilization of engines for subsonic cruise further 
increases the footprint size and enables the vehicle to land as 
a conventional aircraft. Therefore, the amount of crossrange 
and downrange required for a specific mission quite naturally 
depends upon orbit inclination, the location of an acceptable 
landing site, the location of the deorbit maneuver and subsonic 
range capability. In order to define the reentry optimization 
problem, tne following nominal terminal conditions are imposed 

0 r = 6200 mi/r radians 
f e 

= - 1880 mi/r e radians (5.2) 

= 3000 ft/sec 


VI. The Perturbation Method 

The previous sections define the reentry optimization prob- 
lem treated in this study to be of the following form: 


/ 


J 


Q(x,u, t)dt 


( 6 . 1 ) 
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subject to the differential constraints, 


x = f(x,u,t) 


( 6 . 2 ) 


the control variable inequality constraints, 


C(x,u,t) < 0 


(6.3) 


the initial conditions 


x ^o } = x o 
*0 = 0 


(6.4) 


and the terminal conditions 


M(x £ ,t f ) = 0 


(6.5) 


where Q is a scalar function, x is an n-vector of state vari- 
ables, u is an i vector of control variables, C is a q 
vector of control variable constraints, x q is an n-vector of 
initial conditions, t is the initial time and M is a p 

vector of terminal conditions. 

2 

The necessary conditions for a minimum are written as: 


x = H^(x,u,t) 




0 = H u (x,u,X,y,t) 


( 6 . 6 ) 

(6.7) 

( 6 . 8 ) 


where X is an n-vector of time dependent Lagrange multipliers 
associated with the states and H is the variational Hamiltonian, 


with 


H = Q 

+ X T f 

T 

+ p C 

Pi - o 

when 

c t < 0 

> o 

when 

C. = 0 
i 


(6.9) 


0 


( 6 . 10 ) 
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where y is an ^-vector of multipliers associated with the 
control variable inequality constraints and where the matrix 
H uu must be positive definite. 

The following boundary conditions must be satisfied: 

x = x Q at t = t Q (6.11) 

and 

H(t f ) = 0 

A(t f ) = vTM x ( x £> t f) at t = t f (6.12) 

M (x f , t £ ) = 0 


where v is a p-vector of constant multipliers associated with 
the specified terminal conditions. 

Elimination of the control u by using Eqs . (6.8) and 

(6.10) reduces the optimization problem to a two-point boundary 
value problem (TPBVP) which is stated in the following way: 

Find the unknown elements of z and the final time t £ 
which yield the solution of 


z - F (z , t) 


(6.13) 


where 


and 



(6.14) 


(6.15) 


Such that the (n+1) vector, h(z£, t^) , of terminal con- 
straints, composed of Eq. (6.12) vanishes. 

Numerous techniques are available for solution of two-point 
boundary value problems. The technique used in this study is 
the method of perturbation functions (MPF) discussed in Ref. 3 
and 7 . 

To initiate the method, the unknown elements of the vector 
z 0 (in this study the Lagrange multipliers, A ) and the final 
time are selected in some way. A nominal trajectory is then 
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obtained by integrating Eqs . (6.15) from t Q to t £ . Since 

this trajectory seldom satisfies the terminal conditions, a 
correction procedure is required to determine the necessary 
changes in and t £ (6 A q ancI <5t £ ) . This is accomplished 

by considering small perturbations about the nominal path 
through linearization of the nonlinear differential equations 
and terminal conditions as follows: 


where 


6z = A6z 



(2n x 2n matrix) 


(6.16) 


(6 . I 7 ) 


Also, the expected change in the terminal dissatisfaction is 
written as 


where 



(6.18) 


(6.19) 


If it is desired to satisfy the end conditions in one iteration, 
Ah is chosen to be 


Ah = - h 


so that 



( 6 . 20 ) 

( 6 . 21 ) 


The changes 6z f are related to the changes 6z in this line- 
arized system by the fundamental matrix, < K t f> t 0 ) , according 
to 


<Sz £ = 0*(t £ ,t o )6z o 


( 6 . 22 ) 


where obeys the following linear matrix differential equation 

and initial condition 
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= A$(t,t 0 ) 


(2n x 2n matrix) 


(6.23) 


$(t o’V = 


}_ L o_ 
0 ! 1 


(6.24) 


where I is an n x n identity matrix and 0 is an n x n null 
matrix . 

Since all initial states are specified, only the effect of 
perturbations in the initial multipliers must be computed. 
Therefore, only the right half or the last n columns of the fun- 
damental matrix must be integrated. Following this idea, let 


then 


and 


♦ = [*i ■ * 2 ] 

4> ? = A4>2 (2nxn matrix) 

0 
I 

6z f = $ 2 ^f ’ t o-* 6X o 


wv 


(6.25) 


(6.26) 


(6.27) 


By using this result, Eq. (6.21) is rewritten as 



- — 


— 1 

h = 

ah 
3 z 


• 

h 


— — 

f 

- - 


( 6 . 28 ) 


Solution of this system of n + 1 simultaneous linear equa- 
tions yields the changes 6 A q and 6t£ which, if the TPBVP were 
linear, would produce a new vector A q and final time t£ ca- 
pable of satisfying the terminal constraints. However, the non- 
linear nature of most optimization problems requires an iterative 
procedure in which only a portion of the predicted correction is 
applied on each iterate in an effort to maintain the validity of 

the linearization process. The technique used in this study for 

3 

scaling the corrections is the same as that used by Williamson, 

7 

and Lastman and Tapley in which the scale factor is computed 
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such that the magnitude of the correction vector never exceeds 
a prescribed fraction (e.g., 0.30) of the magnitude of the ini- 
tial vector of multipliers, X 

r o 

With the preceding discussion in mind, the basic algorithm 
of the MPF is presented as follows: 

1. Guess nominal starting values for X q and t^ . 

2. Integrate the differential equations for the states 
and the multipliers simultaneously with the linear 
perturbation equations from t Q to t^ . 

3. Evaluate the dissatisfaction in the terminal conditions 
and compute the changes 6 X q and Stj . 

4. Add the scaled changes to the previous values for X 
and t^ and repeat steps 2 through 4 until the dis- 
satisfaction in the terminal conditions is considered 
small . 

VII. Numerical Integration 

The integration of the differential equations for this in- 
vestigation was performed using the variable stepsize Runge- 
Kutta 7(8) formulation developed and described by Fehlberg^. A 
relative single-step truncation error analysis, based on the 
leading term of the truncation error for the seventh-order formu- 
lation in which both the linear and nonlinear equations were 
considered, was used in computing the stepsize. 

The units used in the integration of the differential equa- 
tions are miles, radians and miles per second. The use of these 
units affords a form of normalization of the variables and mul- 
tipliers which tends to aid in the convergence characteristics 
of the problem. All numerical computation was performed in 
single precision on the CDC 6600/6400 computer system at the 
University of Texas at Austin. 

VIII. Application of the MPF and Numerical Results 

The MPF requires the optimization problem to be reduced to 
a TPBVP . Therefore, begin by writing the variational Hamiltonian 
for the reentry as 
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H = 


Q + A V sin Y + A fl - c ?- s ^ cos » + A V cos y sin j 
x r ' 0 r cos <J> <J> r 

+ A v( 


k sin y _ 1 v2 
'V*' r 2 2 pv J * V 'D 


6 r cos (p 

- i pv 2 s*c n ) + 


J i r k cos y V cos y , 1 , ro „ ^ 

+ A y (~ ~YP + — 7 + 2 P vs * C l cos 0 


+ w. V cos y cos H an <}> 1 ,, sin 0, 

P r '2 pVb * L L cos Y J 


+ u fa 

max 


a) (a . - a) + M„(f? - • 

J mm J p max u '• min 


- P) 


where S* = 


S 

m 


(8.1) 


The differential equations for the multipliers associated 
with the states, 



then become 


« 

X 


r 


Q + 

x r 


V cos y cos iff 
r 2 cos <p 


, V cos y sin \p 

% P 


, r 2k sin y 
A y( p 


1 

2 


p r v 2 s*c D ) 


X. ( 
Y 


2k cos y V cos Y . 1 n ^ 

— Vrl — r 2 7 p r VS * C L cos P) 


, r V cos Y COS if tan <fr 1 „ wc , ^ sin p. 
Ail " O q P V U4ut J 

V r z 2 r * L cos 


( 8 . 2 ) 


(8.3) 


= 0 


X, = 


x„ = 


, V cos y cos if sin 4> . V cos y cos if 

0 r cos 2 <f <f r cos 2 <f 


'V 


- Q v " A r sin y - A 


COS Y COS Ip 


cos y sin if 


\\ r ^ -*- 11 r /i q , 'v . 

v r d r cos <f <P 

r -> rk C °S Y . COS Y . 1 o o 

a v pvs*c d - x^c — - -Y r2 + — + j ps*c L cos e) 


+ 


cos Y COS » tan » + 1 s c sin_£ 
V r 2 K * L cos y^ 


(8.4) 

(8.5) 


( 8 . 6 ) 
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, v rn „ v + 1 V sin y cos ip . V sin y sin \p 
- A r V cos y + a Q r“ cos + % r 

^ k cos y ^ r k sin y V sin y-^ 

A V r 2 M vr 2 “ r J 

_ ( V sin y £°S_$_tan_j . 1 vs c sin e.tan Y) 
\\) K r 2 K * L cos y 


(8.7) 


where 


A , = 


and 


A V cos y sin rp _ , V cos y cos ip _ , 
0 r cos <p A <p r 


5 ■ A c% X, s ih y i 


i = l 
4 

Q r “ x c £ [ ( % g ri + % r *i )s i y i ] 

4 

Q v “ A c £ [t%8v i * V i)S l y i ] 


* n r . 5 u 3 . 1 5 

% T = °- 5 k c p p r V 


q. = 3.15 kp'V ’ 15 


V cos y sin tan 4> 
r 

( 8 . 8 ) 

(8.9) 

( 8 . 10 ) 

( 8 . 11 ) 

( 8 . 12 ) 

(8.13) 


P r = - k d p (8.14) 

The partial derivatives of the boundary layer function, , 

are written as 

9 g-i 3gi 9 R n 

g x i = W~ = “95c (8.15) 

Since g. is dependent upon the ratio R /R. > so are i- ts deri- 
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In addition, 
9R 


and 
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9R 


n . Toe v _ • 3 3 j 4 

— = 1.125 k R p p^V q n 
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The optimal control is determined by requiring H a = 0 , 

H 0 = 0 and the matrix H „ given by 
p uu 
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uu 
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aa 
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8 a 
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to be positive definite. 

From these conditions, the optimal roll angle is given by 


sin 6 = 


cos 8 = 


± 


(Aj + A 2 cos 2 Y ) 
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- cos y 
2 


(Aj + A 2 cos 2 y) 


The optimal angle of attack is given by 

1 1 
a = 2 [h + arcs in (y sin n)] 


where 


sin n = (A^ sin 8 - A^ cos y cos 8)/P 


= ( 


^ c^o 


(8.19) 


- 77 < n £ Tf (8.20) 
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and 


2 * c q 0 f 

P = [cos 2 Y (t L S.g. 

jpVS^C^ i=l 


a v vc d 


g . C . 
1 6 1 1 


k) 

c T J 


+ (A cos Y cos 3 - A^ sin 3) 2 ] 


1/2 


( 8 . 21 ) 


The boundary conditions for the TPBVP are summarized as 
follows 


at t = t 

o 

r - r = 400,000 feet 
o e 

0 O = 0 

*0 ‘ 0 

V o = 25,975 feet/second 

Y Q = - 1.5 degrees 

= 0 


at t = t f 

A = 0 (r f free) 

r f 1 

0£ = 6200/r e radians 

(J>£ = - 1880/r e radians 

V£ = 3000 feet/second 

A = 0 (Y f free) 

Tf t 

= 0 C^ £ free) 

H£ = 0 (t£ free) 


( 8 . 22 ) 


(8.23) 


The values for the unknown initial multipliers, A q , 
and the final time, t£ for the trajectory presented are listed 
below 
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A = - 1.54294713 
r o 

A q = - 3.24144967 
0 o 

A, = 3.08551281 x 
*o 

Ay = - 1.52864923 
o 

A = 1.04803561 x 
Y o 

A , = 59.9992066 x 

r o 

t£ = 2039.729 seconds 

• 

Graphical representation of the resultant trajectory is pre- 
sented in Figures 7 through 11. 

The optimal angle of attack and roll angle are plotted 
versus time in Figure 7. Whereas, the maximum roll angle con- 
straint is encountered twice, the angle of attack encounters 
its minimum boundary only once, near the end of the trajectory. 
For the major portion of the trajectory the optimal angle of 
attack is very nearly that for • This is not unusual be- 

cause the specified crossrange is near the maximum crossrange 
for the vehicle and the downrange is considerable, which nor- 
mally requires high values for E . 

The time histories for altitude, velocity and flight path 
angle are presented in Figure 8; while those for downrange 
crossrange and heading angle are presented in Figure 9. 

Referring back to Figure 7, it is seen that the roll angle is 
at its maximum during the initial reentry phase and the first 
altitude peak, which supports the contention that a roll angle 
constraint would reduce skip altitudes. The terminal flight 
path angle is approximately -22°. Although for operational 
reasons a value nearer zero might be more desirable, the ter- 
minal altitude and velocity should be adequate to allow the 
necessary transition to powered flight. 

The reference heating rate, q Q , the dynamic pressure 
and aerodynamic load factor are plotted versus time in Figure 10. 


x 10 
x 10' 
10" 1 
x 10 1 
10 
10 
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The peak heating rate of 124 BTU/ (FT2 -SEC) which occurs at the 
bottom of the first altitude pull-up, is followed by lesser 
peaks occurring each time the minimum altitude in a pull-up 
maneuver is reached. Dynamic pressure and load factor reach 
their maximum values of approximately 266 psf and 1.42 g's, 
respectively, at the end of the trajectory. The maximum load 
factor is well within projected limits for space-rated per- 
sonnel . 


The temperatures on the fuselage panels considered in this 
study (computed from the following equation 


T. = 
i 


Qi U/4 


e s a 


- 460' 


( ° F) 


(8 


where is the surface emissivity and a is the Stephan- 

Boltzman constant) are plotted versus time in Figure 11. 

The maximum temperatures attained are 


T n = 2249° F 

i MAX 

T ? = 1649° F 

Z MAX 

T.. = 1509° F 

"’MAX 

T, = 1607° F 

4 MAX 


The effect of the growth of the turbulent boundary layer on the 
temperatures is evident in the distinct change in character of 
the heating on panel 4 at approximately t = 1200 seconds. The 
heating to each panel exhibits this effect in turn as the tran- 
sition point moves toward the front of the vehicle. 

The value of the total heat load to which the four panels 
are subjected for this high crossrange trajectory is computed 

7 

to be 3.276 x 10 BTU's. As expected, other trajectories with 
lower crossrange have lower heat loads. In particular, trajec- 
tories computed for crossrange on the order of 1700 miles have 
produced heat loads on the order of 2.5 x 10 BTU's. 
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IX. Guidance 

A three-dimensional optimum reentry trajectory found in 
the earlier part of this study was used as the reference tra- 
jectory for the guidance studies. The dynamical model of 
Section III was used and a first-order guidance algorithm was 
implemented. The first-order algorithm implemented was a 
closed loop version of the lambda matrix control algorithm^’ ^ 
The details of the algorithm are well known and will not be 
repeated here, since use of algorithm for shuttle reentry gui- 
dance produced unsatisfactory performance. 

When state errors were introduced, the lambda matrix con- 
trol guidance scheme produced terminal state errors which were 
three to ten times larger than the errors experienced when no 
guidance at all was used. Also, due to the nature of the per- 
turbation equations, first-order control schemes cannot correct 
state errors in longitude at all. This was felt to be a signi- 
ficant shortcoming. 

Typical guidance simulation results are as follows. Start- 
ing with an initial state error in altitude of + 1/2 mile at 
the initial time, the guidance algorithm missed the desired 
terminal point by over 150 miles while the uncorrected error 
produced a trajectory which missed the desired terminal point 
by about 50 miles. Other first-order simulations resulted in 
even larger terminal state errors. 

The programs developed for the guidance study were de- 
veloped with the idea that the initial choice of guidance algo- 
rithm might well be wrong and that the programming should be 
set up in such a way that the guidance algorithm could be 
easily changed. The programming for the guidance study was 
divided into three parts. 

The first part (Program A) generates the nominal reference 
trajectory and stores it on tape. Program A is a reduced ver- 
sion of the trajectory optimization program developed as the 
primary tool for the reentry optimization study. This program 
segment will remain the same regardless of the guidance algo- 
rithm being tested. 
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The second part (Program B) contains routines for input 
of the reference trajectory, routines for the computation of 
the guidance gains (implementation of the particular guidance 
algorithm under study), and routines for the output onto tape 
of the reference trajectory plus the guidance gains informa- 
tion. When changing guidance algorithms, almost all of the 
required programming changes will be in the gain computation 
segment of Program B. 

The third part (Program C) is very short. The program 
consists of input routines, a numerical integrator, and a 
small section of code which produces control changes from 
prc-calculated gains plus known state errors. Execution time 
for this program is of the order of 10-15 seconds for simula- 
tion of a 1600 second reentry. No state variable error pre- 
diction is used. Instead, the nonlinear equations of state are 
integrated from point to point using the corrected controls 
(control values obtained by adding guidance -produced control 
deviations to the nominal control values). After each inte- 
gration step the new state error is computed and the corrected 
control values for this point are produced. 

It is felt that the philosophy employed in designing the 
three programs for the study of the guidance problem will greatly 
facilitate the determination of guidance schemes which will per- 
form adequately for the reentry guidance task. It is likely 
that an extensive study involving the testing of existent gui- 
dance algorithms will be necessary to find an algorithm which 
will perform adequately for shuttle reentry guidance. It is 
possible that the development of new optimal or sub-optimal 
guidance algorithms may be necessary. In any case, the guidance 
programs developed in this study provide a ready tool for the 
execution of such studies. 

X. Conclusions and Recommendations 

The problem of optimal reentry of a shuttle-type vehicle 
lias been considered. In particular, trajectories which mini- 
mized the heat input to the underside of the fuselage and 
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satisfied requirements of down range, terminal velocity and high 
crossrange have been computed. The earth has been assumed to 
be spherical, non-rotating and to possess an inverse square 
gravitational force field. In addition, the atmosphere has been 
considered to be at rest with respect to the earth and to vary 
exponentially with altitude. Control of the vehicle has been 
affected through variation of the angle of attack and roll angle. 
The aerodynamic coefficients of lift and drag have been considered 
to be independent of Mach number and Reynolds number and were 
obtained from Newtonian flow theory. 

The method of perturbation functions (MPF) , a second order 
technique, has been used to generate the trajectories and, al- 
though troubled at times by the sensitivity of the trajectories 
to changes in initial conditions, has proved to be an effective 
technique for generating families of solutions, once an initial 
trajectory has been obtained. 

The principle areas in which this study warrants extension 
are (1) improved aerodynamic and atmospheric models, (2) im- 
proved methods for generating the corrections to the unknown 
initial conditions and (3) investigation of the use of multi- 
shooting or intermediate matching techniques, as opposed to the 
single-shooting method discussed here, in an effort to reduce 
sensitivity problems. 
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Appendix A: Constants 


This appendix contains all relevant constants used in the 
reentry problem. 


C D = .0786 
o 

C n =2.09 
L 

A c = 1.0 x 10" 7 

k = 1.4076519 x 10 16 ft 3 /sec 2 
k c = 17600 P 0 "°‘ 5 V; 3 * 15 BTU/ ft 2 -sec) 

k d = 4.2 x 10" 5 1/ft 

k = 2 - 5 

8 (1.5) • 3 

k = 5.1 x 10 6 

R (2116 x 5280) 1 • 5 

S = 6084 ft 2 
W = 214,861 lbs 

p Q = 2.7 x 10' 3 slugs/ft 3 

r = 3960 miles 
e 

S 1 = 431 ft 2 
s 2 = 928 ft 2 

s 3 = 1306 ft 2 

s 4 = 1408 ft 2 

R t = 1.94 x 10 4 1/ft 
L 1 

R t = 1.2 x 10 4 1/ft 
L 2 
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R. = 7.0 x 
Lj 3 

3 

10° 

l/ft 


R, = 3.5 x 
L 4 

3 

10 

l/ft 


b 1 = 0.1602 

> 

c i = 

0 . 0781 

b 2 = 0.0505 

9 

II 

CnJ 

u 

0 .156 

b 3 = 0.0419 

y 

C 3 = 

0 . 065 

b 4 = 0.0451 

y 

C 4 

0.039 


= - 704.9 


a 2 

2974.9 

a 3 

4952 . 3667 

a 4 = 

4066 . 7 

a 5 = 

1646 . 4 

V 

CN 

li 

263.0667 


Wr 2 


m = -^4 slugs 



a = 4.761 x 10' 13 BTU/ (sec- ft 2 - (°R) 4 ) 

e = 0.8 
s 
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The following is a list accompanied by a brief description 
of each of the subroutines comprising the program. 

SI 101 '1)1* -- Main program which contains the 
basic iteration philosophy of 
the method. 

START --- Called by SHOPDP, sets up all 
relevant constants and initial 
guesses for the unknown multi- 
pli ers and final time. 

SETUP --- Called by PARGEN, initializes 
the matrix DEP which contains 
the states, multipliers and fun- 
damental matrix. 

| PARCHE -- Called by SHOPDP, computes the 

' fundamental matrix by integra- 

! tion of the perturbation equa- 

1 tions and by numerical differ- 

i ences for comparison. 

! PARGEN -- Called by PARCHE and by SHOPDP. 

Computes the fundamental matrix 
by integration of linear pertur- 
bation equations (MPART = 1) 
or by numerical differences 
(MPART = 2) . 

TNTGRT -- Called by PARGEN and by CORVEC. 

Performs integration of equations 
by RK 7(8) (METH = 1) or by 
RK 4 (5) (METH = 2) . 

0UTPT Prints states, multipliers and 

time. Called by RK 7(8) or RK 4(5) 
according to the print increment 
specified by IP in call to INTGRT. 

Converts altitude from miles to 
feet, longitude to downrange in 
miles, latitude to crossrange in 
miles, and velocity from miles per 
second to feet per second. 

DERZ Derivative routine for states and 

multipliers. Called bv RK 7(8) 
or RK 4(5) . 

DERZST -- Derivative routine for states, 

multipliers and linear perturba- 
tion equations. Called by RK 7(8) 
or RK 4 (5) . 

TERMC Evaluates terminal dissatisfaction 

vector and its norm, ROLD. Called 
by SHOPDP. 


DELIC 


LSSDP -- 
C0RVEC - 
ESTPCT - 

DENSIT - 

CUBERT - 
BNDRYL - 
REG 

UNB0UND - 
KH'H I CH - 

FLAGSE - 

PRFIND - 
U0PT --- 


- Sets up the linear system to be 

solved for the required changes 

SA and 6t_p . Called by SHOPDP. 
of 3 

- Called by DELIC. Solves for <5A 0 
and 6 t^ by Gaussian elimina- 
tion. r The matrix input to 
LSSDP is destroyed. 

- Called by DELIC. Scales the cor- 
rections computed in DELIC in a 
method specified by the parameter 
KCOR. See listing. 

- Function called by CORVEC which 
solves for the percent correction 
corresponding to the minimum of a 
parabola fit through three suc- 
cessive values of the norm of the 
terminal dissatisfaction. 

- Evaluates density and its first and 
second derivatives with respect to 
altitude. Exponential atmosphere. 
Called by EQMAT . 

- Function called by UBOUND which 
evaluates the 1/3 power of a func- 
tion . 

- Called by PRFIND and TESTS. Eval- 
uates the boundary layer functions 
and their derivatives. 

- Function which evaluates the nor- 
malized miss distance to the en- 
trance and exit points on a control 
boundary. Called by TESTS. 

- Computes angle of attack for accele- 
ration constraint. Called by UOPT 
but is not active in this deck. 

- Function called by UOPT which sets 
a flag which indicates whether the 
equations of motion are to be eval- 
uated on a constraint boundary or 
not . 

- Sets all flags for integration and 
boundaries. Called by PARGEN, SETUP, 
CORVEC. 

- Evaluates performance index and its 
derivatives. Called by EQMAT. 

- Evaluates the optimal control or 
bounded control and the relevant 
derivatives of the control. Called 
by EQMAT. 
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TESTS Perforins tests to determine if a 

constraint boundary is exceeded. 

Contains iteration philosophy to 
hit entrance and exit point to 
boundary. Called by RK 7(8) or 
RK 4(5) . 

EQMAT --- Called by DERZ or DERZST. Eval- 
uates derivatives of the states 
and multipliers and the non zero 
elements of the matrix of partial 
derivatives for the linear pertur- 
bation equations. 

RK 7(8) - Called by INTGRT. Variable step 

seventh-order Runge-Kutta integra- 
tor using Fehlberg coefficients. 

RK 4(5) - Called by INGRT. Variable step 

fourth-order Runge-Kutta integra- 
tor using Fehlberg coefficients. 

RKC0N Called by START. Sets up coeffi- 

cients for RK 7(8) and RK 4(5). 

Description of Input 

Data is input to the program through the following name- 
lists. 


INTGRT 

AER0 

PIC0N 

IBC 

FBC 

The following is a description of the variables input 
through the namelists. 

INTGRT: 


SMALL Stopping condition for iteration. If 

the norm of the dissatisfaction vector 
is less than SM/ T L, the program terminates. 

PCTN Specified fraction of correction to be 

taken, see KCOR. 

IPQ Print increment for integration of state. 


multipliers and linear perturbation equa- 
tions. If IPQ > 0 stepsize is based 
only on the state and multipliers. If 
IPQ < 0 , stepsize is based on all 
variables . 
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IPR 

T0L 

KC0R 

CL0SE --- 
KPAR 


DELT 

MPART --- 

METH 

ITERS --- 
SIGDIG -- 


ERPS 

DVID 

KBUG 


Print increment when integrating 
state and multipliers only. 

Error tolerance for integration. 

Specifies type of correction scaling 
procedure. See listing of subroutine 
DELIC . 

Tolerance to which entrance and exit 
of control boundaries are satisfied. 


KPAR < 0 Analytic and numerical 
partials are computed, 
program terminates . 

KPAR > 0 Same as above except 
program continues 
iteration . 

KPAR = 0 No action, program con- 
tinues . 


Initial integration stepsize, seconds. 

MPART = 1 Analytic partials are used. 

MPART = 2 Numerical partials are used. 

METH = 1 Integration by RK 7(8). 

METH = 2 Integration by RK 4(5). 

Maximum number of iterations allowed. 

If ICOUNT = ITERS, program stops. 


If the absolute value of relative 
terminal dissatisfaction of a state 
is less than SIGDIG the dissatisfac- 
tion is weighted according to that 
element of WTF. 


If the absolute value of a variable is 
less than ERPS, absolute truncation 
error is used on that variable in com- 
puting an integration stepsize. 

Factor which reduces the predicted 
stepsize to avoid numerous rejection 
of steps. 

KBUG f 0 All namelists are printed. 
KBUG =0 No namelists are printed. 


AER0 : 

CLZER0 -- Coefficient in Newtonian lift equation. 
CDZER0 -- Coefficient in Newtonian drag equation. 

CDL Coefficient in Newtonian drag equation. 

S Vehicle reference area, ft 2 . 
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WT Vehicle weight, lbs. 

AMAX Maximum angle of attack, degrees. 

AMIN Minimum angle of attack, degrees. 

BMX Maximum roll angle, degrees. 

BMN Minimum roll angle, degrees. 

PIC0N : 

A1 Coefficient of performance index. 

A2 Weighting factor for penalty function 

on dynamic pressure. 

A3 Weighting factor for penalty function 

on flight path angle. 

A4 Weighting factor for penalty function 

on reference heating rate. 

AP Panel, areas , ft 2 . 

BJ Coefficients of function of angle of 

attack in heating equation. 

CJ Coefficients of function of angle of 

attack in heating equation. 

RL Transition Reynolds number for panels. 

TMX Temperature used to specify maximum 

reference heating rate allowable. 

QL Maximum acceleration limit (not active). 

CHAPMA -- Chapman heating constant, 17600. 

JTEST JTEST = 0 Iteration to hit control 

boundaries is performed. 

JTEST ^ 0 Boundaries are ignored. 

NPAN Number of panels. 

IBC : 

ZI Initial conditions for states and 

multipliers and final time. 

ZI (1) = Altitude, feet. 

Z I C 2 ) = Downrange, miles. 

ZI(3) = Crossrange, miles. 

ZI(4) = Velocity, feet per second. 

ZI(5) = Flight path angle, degrees. 

ZI(6) = Heading angle, degrees. 

ZI (7) thru 

ZI(12), Multipliers. 

ZI(13) = Final time, seconds. 
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KINPT If desired, multipliers and final 

time can be input in octal. For 
KINPT = 0 read in the multipliers 
following Namelist IRC in a 3020 
format . 

GAM In order to avoid the singularity 

in flight path angle, integration 
is terminated if the magnitude of 
the flight path angle exceeds GAM. 
The time at which this occurs is 
taken as a new final time and the 
program continues. 


FBC : 


ZFN Terminal conditions for states and 

multipliers. The units are the 
same as for ZI in IBC. 

KTC A vector which specifies which ele- 


ments of the ZFN vector are to be 
satisfied. For example, if KTC con- 
sists of 7, 8, 3, 4, 11, 12. Then, 
the terminal conditions to be satis- 
fied are 

A - ZFN (7) = 0 

r 

A Q - ZFN (8) = 0 

<p - ZFN (3) = 0 

V - ZFN (4) = 0 

A y - ZFN (11) = 0 

A^ - ZFN (12) = 0 

as well as 

H(t f ) = 0 

for free final time. 


JTFIX --- JTFIX = 1 

Normal Newton-Raphson 
solution for 6 A and 


6t f . 0 

JTFIX = 2 

Fixed final time. 

JTFIX = 3 

Selected elements of the 
KTC vector are satisfied 


in a least square manner 
according to LSQ. 


LSQ A vector which specifies which of the 

elements of the KTC vector are to be 
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satisfied in a least square manner. 
For example, if LSQ consists of 1, 
2, 5, 6, 7 then for those elements 
specified in KTC above, only termi- 
nal conditions on the multipliers 
* r , X Qf A , and H(t f ) are 

treated . 


NTC Number of terminal conditions speci- 

fied for the least square solution. 
In the example above, NTC = 5 . 

WTF Weighting values to be applied to a 

terminal dissatisfaction to reduce 
its importance in the solution for 
the changes 6X q and 6t^ . 

N0J0Y Print flag. 


|N0J0Y| = 1 Angle of attack and roll 
angle are printed on each 
integrated step. 

N0J0Y < 0 Every iterate to enter or 

exit a control boundary is 
printed . 



